Magnetization profile and core level spectroscopy in a multiply quantized vortex of 

imbalanced Fermi superfluids 
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The core structure of multiply quantized vortices is theoretically investigated in fermionic super- 
fluid near Feshbach resonance. Under population imbalance in two hyperfine spin states, the vortex 
core is filled in by the "paramagnetic moment" . Here, we find the spatial oscillation of the magne- 
tization inside the core sensitively due to the topological structure of the pairing field, in the range 
from the weak coupling regime to the unitary limit. This magnetization inside the giant core reveals 
the winding number of the vortex and directly results from the low-lying quasiparticle states bound 
inside the core. It is therefore proposed that the density profile experiment using phase contrast 
imaging can provide the spectroscopy of novel core level structures in giant vortices. To help the 
understanding on these outcomes, we also derive the analytic solution for the low-lying quasiparticle 
states inside the core of a multiply quantized vortex. 

PACS numbers: 05.30.Fk, 03.75.Hh, 03.75.Ss, 47.32.-y 
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I. INTRODUCTION 

Much attention has been focused on quantized vortices 
of both bosonic and fermionic superfluids. In the bosonic 
system, the order parameter is directly associated with 
the particle density, which is observable via the absorp- 
tion image in experiments. Hence, the vortex core, de- 
fined as the zeros of the order parameter, can be visual- 
ized via the absorption image in the Bose-Einstein con- 
densate (BEC) of ultracold atoms, which has provided an 
opportunity to investigate both the static and dynamic 
properties of vortices [l], 0| • 

In the fermionic case, however, the situation is dras- 
tically changed, where the order parameter corresponds 
to the wave function of the Cooper pair and thus the vi- 
sualization of vortices is not trivial. The particle density 
is sensitively affected by the quasiparticle structure with 
the eigenenergy close to the Fermi energy Ep. The pio- 
neering work in 1964 [U revealed the fact that for a singly 
quantized vortex state in an s-wave case, there exists the 
novel quasiparticle state tightly bound in the core re- 
gion ~£, called the Caroli-de Gennes-Matricon (CdGM) 
state. The important points to note are as follows; (i) 
The eigenenergy of the CdGM state of the vortex core 
is embedded in the low energy region, El Ep = 2/(fcp£) 2 
with kp being the Fermi wave number [3|, and its lowest 
energy level is slightly shifted from the Fermi level. Also, 
(ii) the eigenenergy inside the core is well discretized. 
Hence, the spectrum of the quasiparticle states yields 
the particle-hole asymmetry inside the core 0]. It has 
been found that due to the asymmetric and discretized 
spectrum of the CdGM state, the particle density can 
be suppressed inside the core [5|, called the "quantum 
depletion" , which makes a vortex visible via the density 
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profile even in the Fermi system. The depletion can be 
gradually enhanced as the core radius £ approaches the 
mean inter-particle distance kp 1 , i.e., the quantum limit 
U B S H S El EL El • This fact is directly observed in 
ultracold atomic systems under Feshbach resonance [H| . 

The microscopic studies on multiply quantized or giant 
vortices have been carried out in the type-II superconduc- 
tors [II [II [ll, EL EDI, which indicate that the number of 
the branch of the CdGM state is closely associated with 
the topological structure of the order parameter . Re- 
cently, a macroscopic manifestation of giant vortices in 
neutral fermionic atoms was theoretically studied by Hu 
and Liu [24J, who found the oscillating pattern of the 
particle density profiles inside the core in the population 
balanced system. 

The giant vortex state in ultracold Fermi atoms has 
not been experimentally accomplished so far. In BEC's, 
in contrast, the giant vortices have been experimen- 
tally created by using several techniques: The topolog- 
ical phase imprinting method [l^, [20, E[ , a fast rotat- 
ing BEC confined in a quadratic plus quartic trap [2^ |. 
and the transfer of the orbital angular momentum from 
Laguerre-Gaussian photons to a BEC [23(. In addition, 
the Feshbach resonance, which controls the s-wave scat- 
tering length a, enables one to continuously connect the 
BEC with the Fermionic superfluid, i.e., the BEC to 
Bardeen-Cooper-Schrieffer (BCS) crossover. By sweep- 
ing the magnetic field, the BEC of long-lived molecules in 
a > can be transferred into the Cooper pair of fermionic 
atoms with the weak attractive interaction a < 0, across 
the unitary limit a— >±oo. 

In this paper, we propose a manifestation due to the 
topological structure of the order parameter with a mul- 
tiply quantized vortex. Under population imbalance in 
two hyperfine spin states (f- and |-spins), the occupation 
difference of the CdGM state makes the core of a singly 
quantized vortex with the winding number w = 1 be mag- 
netized [HHf!]. Here, in the BCS regime a < 0, by solving 
the Bogoliubov-de Gennes (BdG) equation under popu- 
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lation imbalance, we find that the magnetization profile 
inside the core of the multiply quantized vortex yields 
an oscillation pattern with several peaks located on a 
concentric circle. This oscillation is well understandable 
from the novel quasiparticle embedded in the vicinity of 
Ef- To help this understanding, we also derive the ana- 
lytic solution for the CdGM states in multiply quantized 
vortices, by following the procedure which was proposed 
by Caroli et al. The analysis is extended to the 

vicinity of the unitary limit (a — * oo) of harmonically 
trapped Fermi gases. Here, we discuss that in the BCS- 
BEC crossover, the oscillation pattern in the magnetiza- 
tion profile inside the core is changed from that in the 
BCS regime due to the strong coupling effect, which em- 
bodies the shift of the energy level of the low- lying CdGM 
state. 

This paper is organized as follows. In Sec. II, af- 
ter introducing the theoretical framework based on the 
BdG equation, we present the numerical results on multi- 
ply quantized vortex states under population imbalance. 
This is carried out in the weak coupling regime in Sec. II. 
The analysis is extended to trapped Fermi gases with 
population imbalance in the vicinity of Feshbach reso- 
nance in Sec. III. Here, we discuss the crossover from 
the weak coupling regime to the Feshbach resonance re- 
gion. Following the previous work [27j , we introduce the 
rcgularization of the gap equation and the particle num- 
ber equation, which allows one to describe the micro- 
scopic structure in the BCS-BEC crossover. Then, we 
shall show the numerical results on the macroscopic and 
quasiparticle structures in resonant Fermi systems with 
giant vortices. The final section is devoted to conclusions 
and discussions. In Appendix, the details on the analyt- 
ical solution of the BdG equation in giant vortices shall 
be described. Throughout this paper, we put 7i=feg = l. 



V(r) denotes a background potential. Here, we impose 
the periodic boundary condition along the z-axis with the 
period L. Since the resulting system yields the axisym- 
metry, the quasiparticle wave function may be expressed 
by 



u q (r) 
v q (r) 



„iqe6 



u q {r)e l 
Vq(r)e 4 (<w- tt ') 



(3) 



where the eigenstate is labeled by the quantum num- 
ber q = (q r ,q e ,q z ) with q g = 0,±1,±2, ••• and q z = 
0,2n/L,4n/L,---. 

The BdG equation is self-consistently coupled with 
the gap equation given by 



A(r)=g ]T u q (r)v*(r)f(E q ), 

\E q \<E c 



(4) 



where f{E q ) = l/{e Eci ' T + 1) is the Fermi distribution 
function, E c is the energy cutoff, and g is the coupling 
constant. Because of the breaking of the particle- hole 
symmetry due to the chemical potential shift Sfi/2 = 
/if — /ij,, the sum in Eq. (U) is carried out for all eigen- 
states with both the positive and negative eigenenergies. 
Throughout this paper, we set Sfi > 0, which means the 
spin-up component becomes majority. 

In this section, in order to focus on the vortex core, 
we restrict our attention to the vortex state with the 
weak pair potential Ao = Q.lEp in the absence of the 
background potential, i.e., V(r) — 0. Here Ao is the 
pair potential in the bulk at T = 0, and Ep is the Fermi 
energy in an ideal Fermi gas. Then, the BdG equation 
can be numerically solved by using the Bessel function 
expansion [3~0 | 



II. GIANT VORTICES WITH POPULATION 
IMBALANCE: WEAK COUPLING LIMIT 

A. Theoretical framework 

Let us consider a single vortex state with arbitrary 
winding number w > 1 , where the pair potential can be 
expressed in the cylindrical coordinate r = (r, 9, z) as 



A(r) = A(r)e i 



(1) 



Without the loss of generality, A(r) is the real function. 
We start with the BdG equation for the quasiparticle 
wave function u a (r) and y a (r) labeled by the quantum 
number q [H M, HE H, 11 , 



/C T (r) A(r) 
A*(r) -JCi(r) 



u q (r) 
v q (r) 



E n 



u q (r) 
v q (r) 



(2) 



The single particle Hamiltonian is given by IC a (r) = 

~Jm ~ ^ a V( r )i where M is the mass of fermions and 
/Ltf^ =/z±(5/Lt is the chemical potential of the spin a state. 



u q {r) 
v q {r) 



N B 

'£ 

i=l 



rr(?e>9*),„(9 



V) 

(Q6,q z ) (qo-w) 



(r) 



,(5) 



where the basis function, 



V2 



RJ„ +1 (a[ u) ) 



Jv 



(") 

21 , 

R 



satisfies 



the 



orthonormal 



Jo" (fif (r)ip^\r)rdr = Sij. Here J v 



(6) 

condition, 
is the v-th 



Bessel function and otf^ is the i-th zero of J u . Through- 
out this section, all physical quantities are scaled by 
the length unit hp 1 = \/2ME~f and the energy unit 
Ep : and we fix the radius of the cylinder R — 200k p 1 , 
the height L = 50k p 1 , the number of the basis function 
7Vb = 100, and the chemical potential fj, = Ep- Note that 
the situation described in this section is applicable not 
only to ultracold Fermi atoms but also superconductors 
under a magnetic field, acting on electron spins. 
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B. Quasiparticle structure and local magnetization 
inside giant cores 

In an isolated vortex with arbitrary winding number 
w > 1, the eigenenergy of the CdGM state, whose wave 
function is localized around the core, is analytically given 
from the BdG equation ([2]) as 

_ / W \ UJ Q ( w — 1 \ . , . 

E " = \ Qe 2 ) Ma) + T + —) Bm( ° )wi ' (?) 

with n — 0, ±1, ±2, • • • and qg — 0, ±1, • • • . We also intro- 
duce sin 2 (a) = 1 — q^/(2M/j). For simplicity, from now 
on, we consider the eigenstates with q z =0, i.e., a — n/2. 
The two coefficients ujo,i are shown in Eq. (|A19|1 and the 
details on the derivation of Eq. j7]) are also described in 
Appendix. 

The expression of E q in Eq. is composed by two 
different energy scales: (i) The energy scale comparable 
with the pair potential, u>\ ~ Ao and (ii) the much smaller 
energy scale wo^ Aq/£f<wi in the weak coupling limit 
Ao <C-E_f- For a single vortex state with w = l, the lowest 
branch of the dispersion relation, i.e., n — 0, reproduces 

A 2 1 

the well-known CdGM eigenstate with E q ~ — (qg — 5) 
The lowest CdGM state yields the "pseudo-zero" 
energy E q ~ A^/Ep <C A in the weak coupling regime. 
The wave function is localized in the core region within 
r<£ (see Fig.QJa)), where £ is the coherence length, £ = 
fcjr/(MAo), and characterizes the core radius of the w= 1 
vortex. In addition, Eq. (JT]) for arbitrary winding number 
W > 1 is in good agreement with the result obtained from 
the semiclassical approximation [i"8| . 

Here, two important facts arising from Eq. (J7|) should 
be mentioned: (i) For the odd number of vorticity w, 
there always exists the core bound CdGM state with the 
pseudo zero energy and qe = 0, because the second term 
in Eq. ([7]) can vanish. In general, the wave function is 
well expressed by Eq. (|A3|) in the region within r <c£, i.e., 
u q (r) ~ Jq B {kpr) ~ (fci?r)' 9e '. Hence, the wave function 
of the lowest CdGM state with qg = may have a large 
intensity on the vortex center. In contrast, the pseudo 
zero "core bound" state with qg — never appears in the 
case of the even number w, because the second term in 
Eq. remains finite. The energy of the low-lying state 
on the vortex center is of the order of Ao. (ii) Apart from 
the core bound state with qg = 0, for w > 1, there may 
exist the state with the pseudo zero energy E q ~ even 
in the even number of w when \qg\^0(kp^). In addition, 
the number of such state is determined by the winding 
number of vortices, as we will show below. Hence, it is 
concluded that the pseudo zero energy states can appear 
inside or around a giant vortex core with arbitrary wind- 
ing number, while the peak position of the wave function 
sensitively depends on the winding number, namely, the 
topological structure of A(r). 

To numerically confirm the direct relation between the 
winding number and the pseudo zero energy states, it 
is worth to visualize the quasiparticle structure of the 



w-fold quantized vortex in balanced systems. Here, we 
introduce the local density of states (LDOS), J\f a (r,E), 
which is given by 

M ] (r,E) = J2\u q (r)\ 2 S(E-E q ), (8a) 

M l (r,E) = J2\v q (r)\ 2 S(E + E q ). (8b) 
1 

Figure 1 shows the LDOS for vortex states with w = 
1, 2, 3, 4 in the case of population balance, 5a — 0, where 
N{r,E) = (Af r (r,E) + Af l {r,E))/2. In the core of the 
odd- number vortex (w = 1,3), the lowest eigenstate is 
situated at r = and in the vicinity of the Fermi level 
(E~0). As we have mentioned above, the vortex state 
with the even w has the distinct energy gap on the cen- 
ter of the vortex, where the excitation spectrum is almost 
symmetric with respect to the Fermi level. The CdGM 
branches in the giant vortex with w > 1 always touch the 
zero energy in the position far from the vortex center, 
e.g., for w — 2, the CdGM branches in Fig. [Ijb) become 
zero around rkp ~ 12, which is comparable with the co- 
herence length /cf£ = 20. In addition, the number of the 
pseudo zero energy state depends on the winding number 
of A(r). For instance, it is seen from Fig. [ljc) that the 
vortex state with w = 3 has two pseudo zero modes: The 
wave function of one mode has a peak at r = and the 
other is at r~25fc^ 1 . Hence, it is numerically confirmed 
that the peak position and the number of the pseudo 
zero energy states depends on the winding number w of 
vortices. 

The presence of the splitting of the Fermi level be- 
tween two spin states, i.e., 5)i>0, enables to reveal the 
two characteristic features of the CdGM states inside the 
giant vortex core, namely, (i) the pseudo-zero or gapful 
excitation at the core r = and (ii) the "pseudo" gapless 
excitation at r~£. Figure [2] shows the local "magnetiza- 
tion" profile, 

m(r) = Pl {r) - Pl {r), (9) 

around the core region of giant vortices with w= 1, 2, 3, 4. 
Here, P|,j.(^) are the local particle densities in hyperfine 
spin states, which are given by 

Pl(r) = J2Mr)\ 2 f(E q ), (10a) 
q 

Pi(r-) = ^K(r)| 2 [l-/(^)]. (10b) 
1 

The magnetization of the vortex core results from the 
fact that the excess atoms of the majority component can 
be accommodated by the eigenstates with the pseudo- 
zero energy. For instance, the core of the single quan- 
tized (w = l) vortex is filled in by the paramagnetic mo- 
ment, as seen in Fig. [Ha). The accommodation of the 
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magnetic moment inside the core is associated with the 
low- lying quasiparticle state, i.e., the CdGM state. Since 

Pu(r) = JZz *\d(r, E)f(E)dE ~ J?*'" M(r, E)dE at 
T = 0, the local magnetization m(r) is determined by the 
eigenstates embedded in the spacing between the Fermi 
surfaces of two spin species. Here, Af(r, E) is the LDOS 
at Sfj, — 0, shown in Fig. [TJ In the case of imbalanced 
population, the Fermi energy of the majority (minority) 
spin component is shifted from that of the balanced case 
upward (downward), Ey) = E F + Sfj (Ef' = E F — Sfj). 
Hence, the low-lying CdGM states are embodied by the 
local paramagnetic moment inside the core, when 5fj = Q. 

For the odd-number vortex, since the lowest CdGM 
state at the vortex center has the much small energy 
~ Aq/Ef , the vortex center is always magnetized in the 
presence of the small splitting of the Fermi level, 5/x<C Ao, 
as seen in Figs. HJJa) and (c). This is contrast to that in 
the even number vortex, where the vortex center can not 
be magnetized in the situation of Sfj Ao because of the 
large energy gap of the CdGM state. Furthermore, it is 
seen from Figs. [UJb)-(d) that the pseudo gapless excita- 
tion at r leads to the concentric oscillation pattern of 
the local magnetization, depending on the winding num- 
ber w. For instance, in the case of the doubly quan- 
tized vortex with w = 2 in Fig. [IJb), since the branches 
of the CdGM state touch the Fermi level at r = l2kp X 
and their excitations have almost zero energy, the excess 
atoms can be accommodated there under small splitting 
of the Fermi level 5/x<C Aq. The peak position of m(r) 
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FIG. 1: (Color online) Local density of states Af(r, E) around 
the core in multiply quantized vortices without population 
imbalance in the case of w = l (a), w — 2 (b), w = 3 (c), and 
w — 4 (d) at T = 0, where only the eigenstates with q z = Q are 
taken into account. The solid line shows the corresponding 
local pair potential ±A(r). Throughout this paper, the Fermi 
surface at <5/x = is put on E = 0. 
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FIG. 2: Local magnetization m(x,y) around the vortex center 
in multiply quantized vortices with w = l (a), 2 (b), 3 (c), and 
4 (d). All the data are at <5/i = 0.02£ F and T=0. 



at r ~ 12kp l in Fig. [2jb) coincides to the spatial profile 
of the pseudo zero states in Fig. [TJb). The w — 4 vor- 
tex has two crossing points of the CdGM branches on 
the Fermi level as seen in Fig. [TJd), leading to the dou- 
ble peak structure on the concentric circles with radius 
rk F ~\5 and 30 in Fig.^d). 

Before turing to the unitary limit, we mention to the 
temperature dependence of the local magnetization in- 
side the giant vortex core. Figure [3] shows the local mag- 
netization around the core region r<50kp l of the w^2> 
vortex at several temperatures, T/Tp — 0,0.01,0.02, 
where Sfj, = Q.02Ep is fixed. It is found that the tem- 




FIG. 3: Local magnetization around the core region of the 
w = 3 vortex within r < 50k' 1 at T/T F = (solid), 0.01 
(dashed), 0.02 (dashed-dotted line). The chemical potential 
shift is fixed as 5fx=0.02E F and max\A(r)\=0.1E F at T=0. 
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perature comparable with S/i smoothes the concentric 
peak structure in m(r). In particular, we should em- 
phasize that the paramagnetic moment in the outside 
of the region r > 30/cp 1 increases as T increases, which 



is characterized by the Yosida function [32[, while that 
at the vortex center (r — 0) and at the concentric peak 
(r~25fc^ 1 ) decreases. 



Ref. [27|. For all data, the total particle number is con- 
served as TV = 9, 000, and we also set L = 3d, correspond- 
ing to the Fermi energy 

Here, the trap length and energy scales are introduced, 
i.e., d = \J 1 /Mui and u), respectively. The calculation is 
carried out in the range from the weak coupling regime 
(fcpa) -1 = — 1.2 to the unitary limit (kpa)^ 1 = 0. 



III. GIANT VORTICES WITH POPULATION 
IMBALANCE: TRAPPED FERMI GASES NEAR 
FESHBACH RESONANCE 

A. BCS-BEC crossover theory 

Having obtained the direct relation between the quasi- 
particle structures and magnetization inside giant vortex 
cores in the weak coupling limit, let us now proceed to ex- 
tend the analysis into the more realistic situation, such as 
a trapped Fermi gas under an s-wave Feshbach resonance. 
Here, we consider the three dimensional cylindrical sys- 
tem that the fermions are confined by the two dimen- 
sional trap potential V(r) = ^Mui 2 r 2 , where r 2 = x 2 + y 2 . 
To do with the vicinity of s-wave Feshbach resonance, 
we have to modify the theoretical framework based on 
the BdG equation in Sec. II A. First, the ultra-violet 
divergence in Eq. ((4|) is removed by replacing g to the ef- 
fective coupling constant g, which is associated with the 
energy cutoff E c and the dimensionless coupling constant 

i/{k F a) mm, 



Ep 
T 

F 



1 



gk 3 



1 



8irkpa 4tt 2 



(11) 



where a is an s-wave scattering length. Also, the chemical 
potential fi is adjusted during the iteration to conserve 
the total particle number, 



N = N ] +N h N U 



(12) 



where the particle density of each spin state is given by 
Eq. (HI. 

The BdG equation ^ is now self-consistently coupled 
with the particle number equation (fl"2"l) and the regu- 
larized gap equation ([4J, where the bare coupling con- 
stant g in Eq. Q is replaced to the effective one g in 
Eq. (TTTj) . The set of equations is free from the energy 
cutoff E c and allows one to describe the qualitative fea- 
ture of T = superfluid phases in the entire range of 
(kpa)- 1 from the BCS ((kpa)- 1 -> -co) to the BEC 
limit {{kp a) -1 — >+oo). 

We numerically solve the gap equation (g]) up to the 
energy 

^(BdG) = i 00oJrj by using the quasiparticle wave 
function obtained from Eq. ([2]). In addition, the higher 
energy contribution above _E( BdG ) is supplemented by 
employing the local density approximation (LDA). The 
details of the BdG-LDA hybrid method are described in 



B. Quasiparticle structure and local magnetization 
inside giant cores 

Here, since in the actual experiment the net magneti- 
zation is conserved due to the absence of the spin relax- 
ation process, it is convenient to introduce the population 
imbalance between two hyperfine spin states, 



N^+N L 



(13) 



Similarly to the non-trapped system in Fig. [21 it is 
found that the distinct pattern of the local "magneti- 
zation" m(r) for vortices with arbitrary winding num- 
ber appears inside the core in the weak coupling regime 
l/(kpa) = —0.4,-1.2, even in the presence of the trap 
potential. Note that the coherence length is estimated 
as f = k F /{MA ) = b^kp 1 = 0.52d at l/(k F a) = -0.4 
and £ = WJkp 1 = 0.167d at l/(k F a) = -1.2, where 
Ao = max |A(r)| at T = 0. The magnetization profiles 
inside the w = l, 2, 3 vortex cores at l/(k F a) = — 0.4 are 
displayed in Figs. Ufa), Hfc), and [He), respectively, where 
the fermions are confined by the harmonic trap. 

As the interaction approaches the unitary limit, how- 
ever, the profile of m(r) is drastically changed. Fig- 
ures [D^b), 0Jd), and Q|f) show the profile of the lo- 
cal magnetization around the giant vortex core at the 
unitary limit, 1/kpa — 0, where the coherence length 
£ = 3.0/c^ 1 = 0.3<i. In the vortex state with w = 1 at 
l/(kpa) = 0, the lowest CdGM state has a large en- 
ergy gap comparable with A , because of A /Ep = 0(1). 
Hence, the energy gap of the lowest CdGM state bounded 
on the vortex center leads to the absence of the magneti- 
zation inside the core even in the high value of P. Note 
that the large population imbalance in the w — 1 vortex in 
Fig. [31(b) results from the magnetization accommodated 
around the edge region of the cloud, r>^. It means that 
the volume of the magnetization m(x,y) inside the core 
within r/d< 1 is almost same in all of Figs. HJa)-(f). 

At l/(kpa) = 0, in the vortex state with w > 1, we 
also find the completely different behavior of m(r) from 
that at 1/(/cf») = —0.4. For instance, in contrast to 
the pattern shown in Figs. d{b) and [He), the core of 
the doubly quantized vortex in Fig. Sfd) is filled in by 
a large amount of the excess fermions, which yields the 
single peak structure around the vortex center. Similarly, 
the small depletion of m(r) appears at r = in the w = 
3 vortex on a resonance, which is contrast to that in 
Figs. (2Kb) and (He). 
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The behavior of m(r) in the core region of various gi- 
ant vortices results from the strong coupling effect with 
A /Ep ~ 0(1) at 1/kpa = 0. To understand this, in 
Fig. 03a), we first plot the energy shift of the peak of 
Af(r = 0, E) at 5^i — 0, as functions of l/(kpa). As seen in 
Fig.QJb), in the weak coupling regime with A /Ep = 0.1 
and fi = E F , the LDOS's A/" u (r = 0,E) of the w = 2 
vortex yield the double peak structures at E ~ ±0.7Ao 
on the vortex center, whose spectral evolution is almost 
symmetric with respect to the Fermi level E = 0. At 
l/kpa = 0, however, the double peaks of A/"-|-,j.(r = 0, E) in 
the w — 2 vortex are shifted upward and their positions 
become asymmetric with respect to the Fermi energy, 
as seen in Fig. Eta) with the solid line. Here the lower 
branch stays nearby the E^O region in the 1/(/cf(i)~0 
region. This shift is understandable from the analytic 
expression of the CdGM state shown in Eq. (JT]), 

= wo ± iwi, (14) 



to m 




■i : i - o 




- 1 I - I 



FIG. 4: Local magnetization m(x, y) around the vortex core 
with the winding number w = 1 (a, b), 2 (c, d), 3 (e, f) at 
T = 0. Figures (a), (c), and (e) are in the BCS side of the 
resonance 1 / (fc^a) = —0.4, and the others are at the resonance 
\/(kpa) =0. Population imbalance is estimated as (a) P — 
0.12 (to = 1), (b) 0.82 (to = 1), (c) 0.014 (to = 2), (d) 0.03 
(w = 2), (e) 0.014 (to = 3), (f) 0.03 (to = 3). 



with ge=0, 9z = 0, and n = 0, — 1. Here, let us recall that 

A 2 

u)q and u>\ are of the order of and A , respectively, 
i.e., loq -C in the weak coupling regime. Then, the 
eigenenergies of two lowest CdGM states in the w = 2 
vortex are symmetric with respect to the Fermi level E = 
0. Indeed, as plotted at \/(kpa) = —1.2 in Fig. EJa), 
the point below (above) E = corresponds to the plus 
(minus) sign of Eq. (|14|) . As approaching the unitary 
limit l/fci?a = 0, however, uq becomes comparable with 
u>\ , which causes the energy shift of the CdGM state with 
— \<jJ\ located inside the Fermi level toward E — 0, e.g., 
E = -0.05A at l/{k F a) = in Fig. gfa). This energy 
shifts of the CdGM states localized on the center of giant 
vortices are also confirmed in Ref. [24[ , where the CdGM 
states shift up across E = as further approaching the 
BEC limit (l/(k F a) -> +oo). 

Figure E£b) shows the LDOS at r = in w — 2 vortex 
state with and without population imbalance at l/kpa = 
0. In a finite population imbalance, since the energy of 
two spin states is shifted downward ("f -spins) or upward 
(|-spins), the lower branch of the two CdGM states is 
occupied (unoccupied) by the majority (minority) spin 
component, which gives rise to the magnetization of the 
vortex center even in the even number winding vortex. 

Figure E( a) also tells us that as l/(kpa) approaches 
the unitarity, the pseudo zero eigenstate, which appears 
in the w = 3 vortex at l/{kpa) — —1.2, is quickly shifted 
upward, e.g., E — 2Quj = 0.57Ao at the unitarity, as seen 
in Fig. [5ja) andEt c )- It indicates that at l/kpa = 0, the 
distinct energy gap may emerge in the vicinity of E = 
in the case of the odd number w. This energy gap leads 
to the suppression of m(r) at r = as seen in Fig. [U4). 
This is contrast to the case of the weak coupling regime. 
For instance, as seen in Fig. Efle), the magnetization in 
the w = 3 vortex at l/{kpa) — —0.4 has a distinct peak 
at r = and also the other peak on the concentric circles. 
This oscillation pattern of m(r) becomes clear as further 
approaching the weak coupling limit (l/(kpd)— ►— oo). 

In addition, we mention that due to the presence of 
the trap potential, finite population imbalance occurs the 
depairing in the vicinity of the edge, in addition to the 
inner region of the core. The numerical results for the 
pairing field and local magnetization shown in Fig. E] re- 
veals the phase separated state between the superfluid 
and spin polarized normal domains at l/fcpa = 0. Simi- 
larly to the non- vortex [27| and the case of singly quan- 
tized vortex [26(, the pairing field in l/kpa<0 yields the 
oscillation in the surface region, that is, the Fulde-Ferrell- 
Larkin-Ovchinnikov-like oscillation. It is found that the 
oscillation of A(r) is proper to imbalanced systems both 
with and without a vortex line in the extensive range 
l/kpa<0. As the interaction touches the unitary limit, 
however, the coexistence area becomes narrow, since the 
oscillation period is reduced to the interparticle distance 
^kp 1 [13. 
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FIG. 5: (a) Eigenenergies of the CdGM states having the 
intensity at r — as functions of l/(kpa). The solid and 
dashed lines denote the CdGM branch in the w = 2 and 3 
vortex states, respectively. The eigenenergies in (a) are scaled 
by Ao, which is given by Ao = max|A(r)|. LDOS at r — 
0, N^i(r, E), of w = 2 (b) and 3 (c) vortex state with and 
without population imbalance. The all dates in (b) and (c) 
are at l/A)Fffl = 0. 



IV. CONCLUDING REMARKS 

In this paper, we have investigated the core structure 
of multiply quantized vortex in imbalanced Fermi sys- 
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FIG. 6: Profile of the local pair potential A(r) (a) and magne- 
tization m(r) (b) at the unitary limit, l/kpa = 0. Population 
imbalance is estimated in vortex states with various winding 
number as P = 0.143 (w = 2, solid line), 0.126 (w — 3, dashed 
line), 0.123 (w = 4, dashed-dotted line). 



terns with and without a trap potential. In conclusion, 
it is found that in the weak coupling regime, as shown 
in Fig. [2l the local magnetization inside the core of mul- 
tiply quantized vortices yields the concentric oscillation 
pattern, which reveals the quasiparticle structure bound 
inside the core, called the CdGM states. The "pseudo- 
zero" energy of the CdGM states in vortices with an odd- 
winding number enables the accommodation of the mag- 
netization on the vortex center. In contrast, the magnetic 
moment is excluded from the vortex center in the case of 
an even winding number. This is because the CdGM 
state has a distinct energy gap of the order of Ao as 
shown analytically in Eq. ([7]) and numerically in Fig. [TJ 

This difference is alternatively understandable from 
the 7r-phase shift of the pairing field [29|,|3J,|35|,|3^,|37|. In 
general, in the presence of the kink structure of the pair 
potential, quasiparticles across the domain wall feel the 
sign change of the pair potential, giving rise to the eigen 
modes with the zero energy, called the Andreev bound 
state [11] or the mid-gap state [13, [H, HI] . The modes 
are strongly bounded around the domain wall. In the 
case of giant vortices with an odd winding number, since 
quasiparticles tracing the path across the vortex center 
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r = always experience the 7r-phase shift of A(r), the 
bound state, called the CdGM state with the pseudo-zero 
energy in the text, appears in the vicinity of the Fermi 
level, whose energy scale is much smaller than the en- 
ergy gap in the bulk region. In contrast, the pseudo-zero 
modes never appear in case of an even winding number, 
where the pair potential dose not change its sign along 
the path across the vortex center. The relation between 
the topological structure of A(r) and the quasiparticle 
state has been shown analytically in Eq. ([7]) and numeri- 
cally in Fig. [1] The pseudo-zero CdGM states can easily 
participate the magnetization, when the Fermi level of 
two spin states is mismatched. The peak structure of the 
local magnetization appears at r — in Fig. [2a) and (c). 
In the even number of winding, however, the distinct en- 
ergy gap arising from the topological reason prevents the 
accommodation of the excess spins at r = 0, as seen in 
Fig. CP) and (d). 

In contrast to the weak coupling regime, where the 
magnetization inside vortex cores can be closely associ- 
ated with the topological structure of the pairing field, it 
has been found that the magnetization profile inside the 
core of various giant vortices is drastically changed in the 
vicinity of Feshbach resonance. The key factor is that two 
energy scales in Eq. ([7]) , ujq and u>i , becomes comparable 
with each other as l/kpa approaches the unitary limit. 
The upward shit of all core bound states gives rise to the 
drastic change of the magnetization profile inside cores 
from that in the BCS limit. 

Finally, we should emphasize that the oscillating pat- 
tern of the magnetization inside the giant core provides 
the spectroscopy of the discretized core structures. As 
seen in Fig. [fj the magnetization profile inside the core 
is drastically changed as varying l/(kpa). Namely, it is 
proposed that the density experiments using phase con- 
trast imaging reveal the spectrum of the core localized 
CdGM states via the magnetization. 
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APPENDIX A: ANALYTIC SOLUTION OF 
CDGM STATES IN GIANT VORTICES 



In this Appendix, we describe the details on the deriva- 
tion of the analytic expression of the CdGM states in 
Eq. (J7J , from the BdG equation (03) without the chemical 
potential shift, <5/i = 0, and trap potential V(r) = 0. Us- 
ing Eqs. (TTJ) and ^ and assuming q z = k^\/l — sin 2 (a) <C 



fc M = ^j2M\(i\ 2 , the BdG equation (OJ) can be rewritten 



in the cylindrical coordinate to 

'w(q g 



[C m T - 2iM Af 2 ] u q = f s 



2ME n 



u q {Al) 



Here we set A = A(r), u q = u q (r)~[u q (r),v q (r)] T , and 

d 2 W ™2 



Id in ,-■>.■■).. 



m 



with m = 



(A2) 

Also, we use the Pauli matri- 



- wq e + — 

ces Ti ; 2,3 and 2x2 unit matrix fo = diag(l, 1). Throughout 
this Appendix, we consider the BCS regime with [i>0. 

Following the procedure proposed by Caroli et al. 0], 
we introduce a radius r c that A(r) = for r <r c . Then, 
the BdG equation (|A1[) can be analytically solved if one 
assume the following conditions: (i) \qg\ -C r c fcp <C &f£, 
(ii) E q <^A , and (iii) E q < |^| 2 sin 2 (a). 

The wavefunction in Eq. (|A1[) is obtained in the range 
r<r r as 



A u J qe (k + (a)r) 
A v J qf) - w (k-(a)r) 



(A3) 



where A u and A v are the arbitrary constant and we set 

E n 



k± (a) = sin(a) ± 



(A4) 



with v fl (a)=k fl sin(a)/M. 

For r > r c , the wave functions are composed of the 

r(i) 



Hankel function Hm and the slow functions <fi(r) vary- 

(A5) 



ing over the order of £ , 



E 

1=1,2 



Then, Eq. IjAlj) can be reduced to 



A 



To 



dr 



E„ 



2M Vfi r 2 



Vi, (A6) 



and ip2(r) oc ip*(r). Here, we set = v fl (a) and tpi = 
<Pi(r). Under the condition (i)-(iii) described above, the 
right hand side of Eq. (|A6|) can be regarded as the small 
perturbation. We now assume the solution of Eq. (|A6|) 

as 



-iip2(r) 



(A7) 



where tp^ (r) is the solution when the right hand side of 
Eq. (|A6|) is neglected, and ^1,2 (^) is the small correction 

to iff\r), i.e., I^i.a (r 
Bie- x ( r '[l,-i] T with 



to ipf\r), i.e., 2(r)| <1. Hence, one finds ip± (r) 



v M) Jo 



A(r')dr'. 



(A8) 
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Since 1^1,2 Ml Eq. (|A7|) can be also expressed as 

»#i(r) 



(A9) 



Within Eq. (|A"7|) . one can find the solution of Eq. (|A"6|) 



w(q 



E n 



2Mv ll r' 



where we set w M = w M (a) , again. 

To get the solution of the BdG equation (|Alj) , the wave 
functions in two different domains, Eq. (|A3j) for r <r c and 
(IA5p for r>r c , are now matched at r~r c . Because of the 
condition (i), |(70| <Cr c fcp<, making use of the asymptotic 

( 1 21 

forms of J v (z) and Hi ' ' (z) in z^> the wave functions 
for r <r c in Eq. (|A3|) is rewritten as 



/ 2M 



A u cos f fc + r + 



„2 _ 1 
% 4 

2k,r 



2qg 



it (Alia) 



In a same way, one finds the another expression from 
Eqs. (|ATTb|l and l|A"l2|i 



E n 



w(q e 



'-) m - q g + w 



2Mv fl r 



1 



+7 - ( n + - | m. 



(A16) 



where n is the integer. The expressions on ip(r) in 
Eqs. (|A15|) and (j A16|) becomes identical when 7 satis- 
fies 



7 = g I m ~ 



1 



7T 

2"- 



(A17) 



The alternative expressions of ip(r) in Eq. (|A10[) and 
Eq. (|A15|) with Eq. (|A17|) should be identical at r = r c . 
Hence, we finally obtain the eigenvalue of the BdG equa- 
tion ©, 



E„ 



qe 



U! 



2 / sin(a) 



10 — 1 



sin(a)wi,(A18) 



2M ( (qg-w) 
A v cos fc_r + 



2 _ 1 
4 



2fc_r 



2g e - 2w + 1 



(Allb) 



with = ^(a) and fc± = fc±(a). Also, Eq. (|A5|) with 
Eq. jMJ) for r>r c > is 



where 



w 



(Al9a) 



2M 



-xO) 



Sie^+W +B 2 e-" 1+(r) 



(A12) 



with 



J7±(r) = Mv^r + 



,2 _ 1 
4 

2MVfjX 



-7r±V(r). (A13) 



To match two expressions of u q (r) in Eqs. (jAllap and 
(|A12[) at r = r c , one should put the coefficients B\ 2 as 



(AU) 



By comparing with Eqs. (jAllap and (|A12|) . one can ob- 
tain the expression of ip as 

. E q w(q e -%) m~qg 
^(r)~-^r+— — — f- + — - — 7T-7. (A15) 
2Mv^r 2 



nk,, 



2M I . 
h) 



-2x(r') dr > 



(A19b) 



To estimate the order of the energy scale of u>o,i, let 
us consider the simplest case of A(r), that is, A(r) = 
Aotanh(r/£). In this situation, one find lo\ ~ §Ao and 

= fff |£ with /" e-^W = e Here, 5 ~ 0(1). 
Hence, the eigenvalue -Eg is composed of two different 
energy scales, such as Ao and Ag/Ep. The expression 
in the case of a singly quantized vortex (w = l) coincides 
to that in Ref. 3], and that with an arbitrary w repro- 
duces the results within the semiclassical approximation 
in Ref. [H|. The further details on Eq. (|A18p is described 
in Sec. II B. We should mention that the integrals in ujq 
and uii depends on the winding number of A(r). 
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